setwd("C:/Users/Mike/Dropbox/To do/Chapter 1/R")

require(ape)
require(phytools)
require(geiger)
require (mvMORPH)
require(devtools)
require(nlme)
require(qpcR)
require(reshape)
require(ouch)
require(calibrate)
require(surface)


tree<- read.nexus("Chapter 1_trees.nexus")
dat<- read.csv("Chapter 1_data_all PCs.csv", row.names=1)
dat2<- read.csv("Chapter 1_data_all PCs.csv")
datl<-list(as.character(dat2[1,1]),as.data.frame(dat2[1,7]))
uncAIC<- vector(mode="numeric", length=0)
results<- vector(mode="numeric", length=0)
##Ecology###
eco <- dat[,7]
names(eco) <- dat$Species


TreeOnly <- setdiff(tree$tip.label,rownames(dat)) #checks to make sure the tips of the phylogeny matches the morpho data

DataOnly <- setdiff(rownames(dat), tree$tip.label) #checks to make sure the tips of the phylogeny matches the morpho data

#########################################
#Prune Tree
########################################

phyloTime <- drop.tip(tree,TreeOnly)

phyloTimeLadderized <- (ladderize(phyloTime)) 
phyloTimeLadderized <- rescale(phyloTimeLadderized, "depth", 1)
plot(phyloTimeLadderized, cex=0.5, no.margin = T) #Plot ladderized and rescaled tree
add.scale.bar() # Add a simple scale bar indicating the scale for the branches in your tree




## make data set up for adding multiple trophic levels
for (i in 1:length(dat2[,7])){
     datl[[i]]<-list(as.character(dat2[i,1]),as.data.frame(dat2[i,7]))
}
## inputing mutliple trophic options - this should just be done once
datl[[14]][[2]]<-data.frame("blue","yellow")
datl[[15]][[2]]<-data.frame("blue","yellow", "green")
datl[[16]][[2]]<-data.frame("blue","yellow", "green")
datl[[25]][[2]]<-data.frame("blue","orange")
datl[[26]][[2]]<-data.frame("blue","yellow", "orange")
datl[[30]][[2]]<-data.frame("blue","green", "orange")
datl[[31]][[2]]<-data.frame("blue","green", "orange")
datl[[37]][[2]]<-data.frame("blue","yellow", "purple")
datl[[38]][[2]]<-data.frame("blue","yellow", "orange", "red")
datl[[57]][[2]]<-data.frame("green","purple")
datl[[85]][[2]]<-data.frame("blue","yellow")
datl[[92]][[2]]<-data.frame("blue","green", "purple")
datl[[93]][[2]]<-data.frame("blue","green", "purple")
datl[[94]][[2]]<-data.frame("blue","green", "purple")
datl[[95]][[2]]<-data.frame("blue","purple")
datl[[100]][[2]]<-data.frame("blue","red")
datl[[103]][[2]]<-data.frame("blue","red")
datl[[104]][[2]]<-data.frame("blue","red")
datl[[105]][[2]]<-data.frame("blue","red")
datl[[106]][[2]]<-data.frame("blue","green", "yellow", "orange")
datl[[107]][[2]]<-data.frame("blue","green", "yellow", "orange")
datl[[113]][[2]]<-data.frame("yellow", "orange", "blue")
datl[[114]][[2]]<-data.frame("yellow", "orange", "blue")



# ACTUALLY DOING IT
for (j in 1:1000){

# first make dataset to run through programs
# pull ONE trophic classification per species, randomly
     for (i in 1:length(datl)){
          uncNum <- length(datl[[i]][[2]])
          r<-sample(1:uncNum,1)
     dat[i,6]<-datl[[i]][[2]][r]
          }
# output is modified dat with trophic info in col 7
# plug this dataset - dat - into your models

####model for trophic ecology####
discTraitOne <- as.vector(dat[,6])  # Make a vector of discreate trait one.
names(discTraitOne) <- rownames(dat) # Name the vector elements.
cols <- setNames(palette()[1:length(unique(discTraitOne))],sort(unique(discTraitOne)))          

mtrees<-make.simmap(phyloTimeLadderized,discTraitOne,model="SYM", nsim=100)



#####RUN MODEL#####
OUtrophic<-mvOU(mtrees[[1]], dat[,1:4], error = NULL, model = c("OUM", "OU1"), param = list(sigma = NULL,
          alpha = NULL, vcv = "mvmorph", decomp = c("diagonal")), method = c("rpf", "sparse", "inverse",
          "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B",
          "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL,
          diagnostic = TRUE, echo = TRUE)






uncAIC[j]<-OUtrophic$AICc


}



 OUcont =-1735
     OUfam= -1838
     OUpiscivore = -1823
     OUsurface = -1910
     OUsurface_pisc = -1901
     BM = -1769
     EB = -1782

for (i in 1:length(uncAIC)) {
     OU trophic <- uncAIC[i]

results[i]<- list(OUtrophic, OUcont, OUfam, OUpiscivore, OUsurface, OUsurface_pisc, BM, EB)
aicw(results, aicc=TRUE)
}

results<- list(OUtrophic, OUcont, OUfam, OUpiscivore, OUsurface, OUsurface_pisc, BM, EB)

aicw(results, aicc=TRUE)